options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -3993.9 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -38.4782    0.00113674    0.00127173           1           1       47    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -38.48
##    b0         8.80
##    b1         1.66
##    s          4.15
##    m0[1]     18.32
##    m0[2]     18.32
##    m0[3]     31.00
##    m0[4]     46.87
##    m0[5]     31.73
##    m0[6]     40.71
##    m0[7]     20.52
##    m0[8]      7.02
##    m0[9]     38.01
##    m0[10]    37.52
##    m0[11]    38.19
##    m0[12]    33.96
##    m0[13]    38.43
##    m0[14]    22.24
##    m0[15]    34.45
##    m0[16]    17.56
##    m0[17]    39.50
##    m0[18]     7.40
##    m0[19]    29.33
##    m0[20]    21.23
##    y0[1]     19.40
##    y0[2]     16.97
##    y0[3]     23.26
##    y0[4]     43.19
##    y0[5]     31.43
##    y0[6]     38.78
##    y0[7]     18.25
##    y0[8]     11.78
##    y0[9]     41.81
##    y0[10]    41.95
##    y0[11]    33.05
##    y0[12]    36.17
##    y0[13]    36.94
##    y0[14]    19.02
##    y0[15]    35.66
##    y0[16]    16.87
##    y0[17]    33.72
##    y0[18]     5.18
##    y0[19]    22.95
##    y0[20]    22.64
##    m1[1]      7.02
##    m1[2]     11.45
##    m1[3]     15.88
##    m1[4]     20.31
##    m1[5]     24.73
##    m1[6]     29.16
##    m1[7]     33.59
##    m1[8]     38.02
##    m1[9]     42.44
##    m1[10]    46.87
##    y1[1]      6.69
##    y1[2]     11.31
##    y1[3]     21.44
##    y1[4]     23.48
##    y1[5]     22.31
##    y1[6]     28.26
##    y1[7]     32.13
##    y1[8]     44.58
##    y1[9]     41.09
##    y1[10]    50.81
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -38.76 -38.43 1.35 1.15 -41.46 -37.25 1.00      482      871
##    b0       8.73   8.74 2.30 2.35   4.72  12.40 1.02      248      306
##    b1       1.67   1.66 0.17 0.16   1.39   1.95 1.02      303      443
##    s        4.77   4.64 0.87 0.79   3.58   6.31 1.00     1275     1063
##    m0[1]   18.26  18.25 1.53 1.53  15.66  20.66 1.01      276      398
##    m0[2]   18.26  18.25 1.53 1.53  15.66  20.66 1.01      276      398
##    m0[3]   30.96  30.96 1.14 1.09  29.07  32.77 1.00     1681     1317
##    m0[4]   46.85  46.94 2.15 2.06  43.18  50.32 1.01      881      925
##    m0[5]   31.69  31.70 1.16 1.12  29.76  33.54 1.00     1935     1315
##    m0[6]   40.68  40.72 1.65 1.58  37.89  43.28 1.00     1514     1213
##    m0[7]   20.46  20.46 1.39 1.40  18.20  22.65 1.01      301      539
##    m0[8]    6.96   6.96 2.46 2.51   2.64  10.88 1.02      248      284
##    m0[9]   37.98  38.00 1.46 1.39  35.53  40.27 1.00     1991     1448
##    m0[10]  37.49  37.51 1.43 1.36  35.10  39.74 1.00     2060     1517
##    m0[11]  38.16  38.18 1.47 1.40  35.69  40.47 1.00     1963     1384
##    m0[12]  33.92  33.95 1.24 1.18  31.89  35.85 1.00     2330     1362
##    m0[13]  38.40  38.42 1.49 1.42  35.89  40.72 1.00     1928     1384
##    m0[14]  22.19  22.18 1.29 1.30  20.06  24.21 1.01      335      668
##    m0[15]  34.42  34.44 1.26 1.20  32.35  36.38 1.00     2351     1343
##    m0[16]  17.50  17.48 1.59 1.60  14.83  20.00 1.01      269      392
##    m0[17]  39.47  39.50 1.56 1.50  36.86  41.90 1.00     1736     1247
##    m0[18]   7.33   7.33 2.42 2.47   3.07  11.19 1.02      248      284
##    m0[19]  29.29  29.30 1.12 1.10  27.45  31.06 1.00     1009     1254
##    m0[20]  21.18  21.17 1.35 1.36  18.95  23.27 1.01      314      600
##    y0[1]   18.28  18.41 5.12 5.03   9.91  26.66 1.00     1752     1891
##    y0[2]   18.21  18.12 5.18 4.92   9.81  26.80 1.00     1460     1587
##    y0[3]   31.01  30.81 4.93 4.67  23.08  38.91 1.00     1662     2013
##    y0[4]   46.94  46.80 5.42 5.13  38.01  55.79 1.00     1913     1445
##    y0[5]   31.66  31.68 4.94 4.78  23.19  39.70 1.00     2037     1892
##    y0[6]   40.89  41.01 5.04 5.00  32.48  49.15 1.00     1822     1636
##    y0[7]   20.49  20.64 5.09 4.77  11.87  28.60 1.00     1675     1574
##    y0[8]    6.75   6.68 5.41 5.36  -2.08  15.80 1.00      902     1620
##    y0[9]   37.97  38.05 5.11 5.01  29.77  45.99 1.00     1864     1972
##    y0[10]  37.51  37.55 4.95 4.61  29.43  45.49 1.00     1945     1695
##    y0[11]  38.34  38.26 5.08 4.75  30.37  46.57 1.00     2066     1918
##    y0[12]  33.95  33.96 5.01 4.80  25.82  42.55 1.00     2124     1889
##    y0[13]  38.48  38.36 5.16 4.85  30.11  47.08 1.00     1968     1788
##    y0[14]  22.20  22.22 4.99 4.60  13.93  30.51 1.00     1823     1725
##    y0[15]  34.23  34.24 5.01 4.86  26.02  42.33 1.00     1990     1947
##    y0[16]  17.82  17.75 5.07 4.88   9.96  26.42 1.00     1589     1666
##    y0[17]  39.45  39.48 5.06 4.79  31.26  47.62 1.00     2089     1780
##    y0[18]   7.03   7.15 5.35 5.08  -1.65  15.63 1.01      689     1354
##    y0[19]  29.17  29.31 5.06 4.73  20.52  37.06 1.00     1941     1730
##    y0[20]  21.24  21.22 5.07 4.66  13.03  29.38 1.00     1814     1893
##    m1[1]    6.96   6.96 2.46 2.51   2.64  10.88 1.02      248      284
##    m1[2]   11.39  11.37 2.07 2.14   7.77  14.65 1.02      250      309
##    m1[3]   15.82  15.81 1.71 1.73  12.92  18.52 1.02      260      358
##    m1[4]   20.25  20.25 1.40 1.40  17.97  22.44 1.01      299      507
##    m1[5]   24.69  24.69 1.19 1.18  22.75  26.54 1.01      428      905
##    m1[6]   29.12  29.13 1.12 1.10  27.29  30.91 1.00      967     1254
##    m1[7]   33.55  33.58 1.22 1.16  31.56  35.48 1.00     2292     1362
##    m1[8]   37.98  38.01 1.46 1.39  35.54  40.28 1.00     1990     1448
##    m1[9]   42.42  42.47 1.78 1.70  39.36  45.29 1.00     1274     1149
##    m1[10]  46.85  46.94 2.15 2.06  43.18  50.32 1.01      881      925
##    y1[1]    6.99   6.94 5.40 5.32  -1.91  15.84 1.00     1071     1128
##    y1[2]   11.22  11.28 5.37 5.05   2.46  19.87 1.00     1063     1647
##    y1[3]   15.87  15.99 5.07 5.07   7.55  24.14 1.00     1427     1553
##    y1[4]   20.38  20.35 4.97 4.82  12.23  28.50 1.00     1529     1869
##    y1[5]   24.42  24.38 5.00 4.80  16.40  32.71 1.00     1634     1865
##    y1[6]   29.09  29.06 4.95 4.68  21.19  37.12 1.00     1910     2057
##    y1[7]   33.65  33.72 5.07 4.89  25.47  41.85 1.00     2055     1999
##    y1[8]   37.92  38.04 5.09 4.95  29.57  46.02 1.00     2092     1931
##    y1[9]   42.41  42.33 5.08 5.07  34.06  50.75 1.00     1869     1790
##    y1[10]  46.84  46.85 5.34 5.38  38.32  55.71 1.00     1915     1856

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.4 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 finished in 0.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -47.63 -49.21 10.77 9.41 -63.07 -26.53 1.07       54      101
##    b0       8.98   8.86  2.37 2.36   5.43  12.85 1.03      143     1209
##    b1       1.65   1.65  0.18 0.18   1.37   1.93 1.04       80     1155
##    s        3.15   3.24  1.54 1.67   0.70   5.62 1.09       32       26
##    sx       1.98   2.00  0.96 0.97   0.48   3.68 1.07       59       70
##    x0[1]    3.66   3.81  1.97 2.06  -0.18   6.37 1.05       73       59
##    x0[2]    7.02   7.03  1.47 1.62   4.73   9.37 1.03      130     2385
##    x0[3]   13.80  13.75  1.23 1.03  11.88  15.88 1.02     1246     1828
##    x0[4]   21.65  21.86  1.61 1.61  19.08  24.01 1.04       69     1274
##    x0[5]   15.01  14.97  1.49 1.53  12.76  17.47 1.02      192     1391
##    x0[6]   20.72  20.54  1.70 1.66  18.29  23.70 1.04       94      417
##    x0[7]    5.70   5.79  1.63 1.60   2.87   8.09 1.04       97      448
##    x0[8]    0.03  -0.13  1.57 1.58  -2.35   2.68 1.04       77      380
##    x0[9]   16.83  16.88  1.31 1.28  14.79  18.91 1.02      168     1896
##    x0[10]  16.60  16.67  1.30 1.28  14.52  18.69 1.03      166     1831
##    x0[11]  18.39  18.22  1.43 1.27  16.30  20.82 1.02      244     1596
##    x0[12]  15.74  15.66  1.26 1.13  13.80  17.87 1.01      862     1525
##    x0[13]  17.24  17.33  1.29 1.21  15.24  19.28 1.02      315     2283
##    x0[14]   5.02   5.14  2.47 2.85   0.50   8.46 1.06       55       55
##    x0[15]  16.42  16.29  1.40 1.31  14.25  18.74 1.02      299     1537
##    x0[16]   4.96   5.09  1.38 1.18   2.64   7.05 1.01      597     1326
##    x0[17]  20.15  19.93  1.77 1.78  17.72  23.32 1.04      108      788
##    x0[18]   0.41   0.26  1.64 1.67  -2.06   3.10 1.04       65      659
##    x0[19]  11.12  11.16  1.52 1.62   8.65  13.38 1.03      130      663
##    x0[20]   7.06   7.13  1.31 1.10   4.84   9.12 1.01      675     1888
##    m0[1]   15.14  14.96  2.97 3.50  10.89  20.09 1.05       61      771
##    m0[2]   20.65  20.87  2.47 2.61  16.49  24.25 1.02      230     2675
##    m0[3]   31.77  31.90  2.03 1.78  28.42  35.03 1.01     1354     2562
##    m0[4]   44.66  44.30  2.56 2.42  40.92  49.20 1.02      336     2502
##    m0[5]   33.75  33.88  2.44 2.66  29.73  37.30 1.02      142     2299
##    m0[6]   43.13  43.35  2.71 3.00  38.61  46.90 1.03       97     2052
##    m0[7]   18.48  18.27  2.53 2.62  14.72  22.73 1.03       93     2255
##    m0[8]    9.16   9.48  2.63 2.53   4.56  12.89 1.01      426     2137
##    m0[9]   36.75  36.56  2.17 1.92  33.36  40.43 1.01      648     2131
##    m0[10]  36.37  36.20  2.09 1.87  33.15  39.84 1.01      600     2420
##    m0[11]  39.32  39.52  2.27 2.02  35.48  42.92 1.02      311     2339
##    m0[12]  34.96  35.11  2.07 1.91  31.40  38.15 1.02     1116     2676
##    m0[13]  37.43  37.30  2.11 1.84  34.09  41.00 1.01     1135     2405
##    m0[14]  17.37  17.28  3.83 5.02  12.06  23.42 1.07       46      329
##    m0[15]  36.06  36.24  2.27 2.26  32.31  39.45 1.01      331     2927
##    m0[16]  17.26  17.27  2.14 1.80  13.81  20.76 1.02     2866     2173
##    m0[17]  42.18  42.37  2.77 3.17  37.59  46.10 1.03       90     2191
##    m0[18]   9.79  10.13  2.74 2.70   5.11  13.58 1.02      278     1915
##    m0[19]  27.37  27.26  2.44 2.57  23.80  31.40 1.03      118     2156
##    m0[20]  20.71  20.61  2.08 1.75  17.26  24.16 1.02     1848     2395
##    y0[1]   15.18  14.23  4.67 4.23   8.90  23.66 1.02      253     2565
##    y0[2]   20.57  21.28  4.30 3.48  12.95  26.79 1.02     1027     2648
##    y0[3]   31.83  32.11  3.96 3.23  25.00  38.20 1.02     3404     2528
##    y0[4]   44.61  43.98  4.26 3.47  38.35  52.04 1.03     1576     2989
##    y0[5]   33.82  34.45  4.27 3.48  26.11  39.94 1.01      951     2539
##    y0[6]   43.12  43.84  4.42 3.59  35.17  49.38 1.01      456     2510
##    y0[7]   18.53  17.94  4.31 3.47  12.11  26.21 1.02      502     2915
##    y0[8]    9.14   9.70  4.34 3.55   1.32  15.72 1.01     1650     2552
##    y0[9]   36.69  36.42  4.07 3.19  30.38  43.95 1.02     3185     2333
##    y0[10]  36.27  35.95  4.19 3.44  29.87  43.39 1.02     3633     2931
##    y0[11]  39.37  39.74  4.21 3.31  31.96  45.87 1.02     2197     2842
##    y0[12]  34.97  35.27  4.01 3.15  28.11  41.29 1.02     2681     2709
##    y0[13]  37.36  37.08  4.12 3.30  30.95  44.33 1.03     2693     2712
##    y0[14]  17.40  16.47  5.14 5.25  10.81  26.89 1.04       84     1993
##    y0[15]  35.99  36.50  4.28 3.44  28.40  42.48 1.02     1241     2159
##    y0[16]  17.33  17.24  4.13 3.24  10.59  24.29 1.03     3836     2231
##    y0[17]  42.11  42.83  4.51 3.91  33.75  48.34 1.01      384     2462
##    y0[18]   9.84  10.45  4.51 3.61   1.75  16.44 1.01     1202     3066
##    y0[19]  27.46  26.83  4.26 3.56  21.19  35.23 1.01      754     2806
##    y0[20]  20.70  20.56  4.11 3.24  14.00  27.50 1.03     3705     3174
##    m1[1]    7.22   7.09  2.54 2.53   3.42  11.33 1.03      131     1268
##    m1[2]   11.62  11.51  2.11 2.11   8.42  15.06 1.03      166     1368
##    m1[3]   16.02  15.93  1.71 1.72  13.42  18.81 1.02      224     1370
##    m1[4]   20.42  20.36  1.37 1.32  18.29  22.63 1.01      493     1576
##    m1[5]   24.82  24.81  1.12 1.04  22.94  26.56 1.02     1269     2009
##    m1[6]   29.22  29.27  1.05 0.95  27.44  30.86 1.01     1824     2137
##    m1[7]   33.62  33.61  1.19 1.20  31.67  35.48 1.01      401     2122
##    m1[8]   38.02  38.00  1.48 1.59  35.61  40.28 1.03      138     2036
##    m1[9]   42.42  42.37  1.85 2.00  39.47  45.26 1.03      106     1685
##    m1[10]  46.82  46.76  2.26 2.41  43.31  50.40 1.04       94      959
##    x1[1]   -1.11  -1.09  2.24 1.67  -4.92   2.69 1.02     4230     2272
##    x1[2]    1.55   1.55  2.17 1.61  -2.05   5.17 1.02     3986     1588
##    x1[3]    4.24   4.24  2.20 1.61   0.59   7.76 1.02     3708     2226
##    x1[4]    7.00   6.99  2.18 1.67   3.50  10.64 1.02     4243     3124
##    x1[5]    9.57   9.60  2.15 1.60   6.01  13.08 1.02     4066     2474
##    x1[6]   12.31  12.26  2.25 1.67   8.75  16.05 1.02     3840     2435
##    x1[7]   14.85  14.87  2.13 1.65  11.30  18.34 1.02     3955     2330
##    x1[8]   17.52  17.56  2.24 1.66  13.70  21.16 1.02     4001     2018
##    x1[9]   20.20  20.20  2.21 1.68  16.67  23.84 1.02     3772     2267
##    x1[10]  22.86  22.87  2.27 1.64  19.18  26.49 1.02     3524     2260
##    y1[1]    7.26   7.10  4.35 4.14   0.31  14.28 1.01      749     2586
##    y1[2]   11.68  11.54  4.06 3.65   5.13  18.27 1.01      944     2129
##    y1[3]   15.95  15.90  3.95 3.29   9.58  22.41 1.01     1321     2520
##    y1[4]   20.38  20.31  3.80 2.94  14.31  26.79 1.02     3401     2765
##    y1[5]   24.78  24.87  3.69 2.85  18.53  30.70 1.03     3594     2582
##    y1[6]   29.24  29.30  3.69 2.83  22.96  35.15 1.03     3755     2313
##    y1[7]   33.56  33.75  3.76 2.81  27.12  39.48 1.02     2704     2455
##    y1[8]   38.08  38.12  3.87 3.17  31.78  44.49 1.01     1790     3036
##    y1[9]   42.34  42.40  3.96 3.44  35.52  48.57 1.01      916     2423
##    y1[10]  46.82  46.96  4.20 3.91  40.05  53.55 1.01      515     2761

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -60267.9 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       26      -33.7153   0.000295205   0.000111983           1           1       71    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -33.72
##    b0         6.00
##    b1         1.86
##    s          3.27
##    m0[1]     11.57
##    m0[2]     21.43
##    m0[3]      8.68
##    m0[4]     19.77
##    m0[5]     13.48
##    m0[6]     13.99
##    m0[7]     17.35
##    m0[8]     13.47
##    m0[9]     19.21
##    m0[10]     7.15
##    m0[11]     9.82
##    m0[12]     6.26
##    m0[13]    17.35
##    m0[14]    15.68
##    m0[15]    20.82
##    m0[16]    13.73
##    m0[17]    20.82
##    m0[18]    17.98
##    m0[19]    19.16
##    m0[20]    20.26
##    y0[1]     12.23
##    y0[2]     24.51
##    y0[3]     10.42
##    y0[4]     24.35
##    y0[5]     12.97
##    y0[6]     18.25
##    y0[7]     18.06
##    y0[8]     12.42
##    y0[9]     20.31
##    y0[10]     7.74
##    y0[11]    10.88
##    y0[12]     1.30
##    y0[13]    15.88
##    y0[14]    18.45
##    y0[15]    24.94
##    y0[16]    22.56
##    y0[17]    19.99
##    y0[18]    19.73
##    y0[19]    19.53
##    y0[20]    18.98
##    m1[1]      6.26
##    m1[2]      7.95
##    m1[3]      9.63
##    m1[4]     11.32
##    m1[5]     13.00
##    m1[6]     14.69
##    m1[7]     16.37
##    m1[8]     18.06
##    m1[9]     19.74
##    m1[10]    21.43
##    y1[1]      5.93
##    y1[2]      8.75
##    y1[3]      6.32
##    y1[4]      9.92
##    y1[5]     16.23
##    y1[6]     13.25
##    y1[7]     12.10
##    y1[8]     19.00
##    y1[9]     15.78
##    y1[10]    19.57
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -34.08 -33.75 1.27 1.05 -36.60 -32.69 1.00      650      951
##    b0       6.10   6.03 1.76 1.76   3.35   8.96 1.01      394      451
##    b1       1.84   1.84 0.32 0.31   1.31   2.36 1.01      474      508
##    s        3.71   3.62 0.68 0.61   2.76   4.96 1.00     1268     1136
##    m0[1]   11.62  11.62 1.02 0.99   9.96  13.33 1.01      456      692
##    m0[2]   21.39  21.38 1.37 1.27  19.11  23.67 1.00     1030     1048
##    m0[3]    8.76   8.73 1.37 1.36   6.63  11.00 1.01      400      472
##    m0[4]   19.75  19.76 1.15 1.03  17.82  21.68 1.00     1292     1152
##    m0[5]   13.51  13.52 0.87 0.83  12.08  14.96 1.01      604     1042
##    m0[6]   14.02  14.02 0.84 0.79  12.62  15.41 1.01      678     1157
##    m0[7]   17.35  17.37 0.91 0.83  15.83  18.84 1.00     1622     1339
##    m0[8]   13.51  13.52 0.87 0.83  12.08  14.96 1.01      603     1042
##    m0[9]   19.19  19.20 1.09 0.97  17.38  21.01 1.00     1412     1147
##    m0[10]   7.24   7.20 1.59 1.58   4.78   9.82 1.01      394      472
##    m0[11]   9.89   9.87 1.22 1.22   7.97  11.87 1.01      412      575
##    m0[12]   6.36   6.31 1.72 1.72   3.67   9.15 1.01      393      438
##    m0[13]  17.35  17.38 0.91 0.83  15.83  18.85 1.00     1622     1339
##    m0[14]  15.70  15.71 0.83 0.77  14.30  17.04 1.00     1099     1304
##    m0[15]  20.79  20.78 1.29 1.17  18.65  22.93 1.00     1116     1151
##    m0[16]  13.76  13.76 0.85 0.81  12.35  15.19 1.01      638     1090
##    m0[17]  20.79  20.79 1.29 1.17  18.65  22.93 1.00     1116     1151
##    m0[18]  17.97  18.00 0.96 0.87  16.36  19.56 1.00     1618     1297
##    m0[19]  19.14  19.15 1.08 0.97  17.34  20.96 1.00     1422     1147
##    m0[20]  20.23  20.24 1.21 1.10  18.19  22.26 1.00     1204     1151
##    y0[1]   11.70  11.68 3.95 3.83   5.30  18.30 1.00     1797     1864
##    y0[2]   21.42  21.34 3.95 3.70  15.14  28.01 1.00     1739     1733
##    y0[3]    8.81   8.75 4.04 3.85   2.13  15.62 1.00     1378     2084
##    y0[4]   19.68  19.71 3.92 3.83  13.26  25.85 1.00     1982     1969
##    y0[5]   13.52  13.63 3.84 3.80   6.99  19.63 1.00     1859     1717
##    y0[6]   14.13  14.07 3.80 3.78   8.18  20.36 1.00     1521     1762
##    y0[7]   17.21  17.26 3.78 3.64  10.67  23.42 1.00     2137     1970
##    y0[8]   13.39  13.52 3.91 3.79   6.80  19.59 1.00     1805     1670
##    y0[9]   19.15  19.11 4.05 3.95  12.66  25.75 1.00     2176     2058
##    y0[10]   7.29   7.31 4.07 3.89   0.64  13.89 1.00     1194     1570
##    y0[11]   9.89   9.88 4.01 3.73   3.22  16.73 1.00     1330     1843
##    y0[12]   6.26   6.21 4.15 4.02  -0.42  12.98 1.00     1220     1416
##    y0[13]  17.16  17.11 3.94 3.81  10.64  23.68 1.00     2149     1975
##    y0[14]  15.70  15.65 3.86 3.65   9.44  22.19 1.00     2003     1949
##    y0[15]  20.86  20.94 4.01 3.91  14.43  27.53 1.00     2020     1917
##    y0[16]  13.81  13.73 4.02 3.94   7.04  20.37 1.00     1913     1362
##    y0[17]  20.76  20.84 4.04 3.82  14.11  27.43 1.00     2064     1892
##    y0[18]  18.02  18.04 3.93 3.89  11.41  24.32 1.00     2005     1770
##    y0[19]  18.99  19.01 4.07 3.86  12.25  25.48 1.00     1987     1955
##    y0[20]  20.23  20.29 3.87 3.75  13.88  26.53 1.00     1735     1857
##    m1[1]    6.36   6.31 1.72 1.72   3.67   9.15 1.01      393      438
##    m1[2]    8.03   8.01 1.48 1.48   5.75  10.42 1.01      397      473
##    m1[3]    9.70   9.68 1.25 1.26   7.75  11.74 1.01      409      569
##    m1[4]   11.37  11.37 1.05 1.02   9.68  13.11 1.01      447      697
##    m1[5]   13.04  13.05 0.90 0.87  11.58  14.56 1.01      550      938
##    m1[6]   14.71  14.72 0.82 0.78  13.36  16.09 1.00      805     1228
##    m1[7]   16.38  16.41 0.85 0.79  14.97  17.77 1.00     1368     1319
##    m1[8]   18.05  18.08 0.97 0.88  16.42  19.65 1.00     1605     1297
##    m1[9]   19.72  19.73 1.15 1.03  17.80  21.65 1.00     1299     1152
##    m1[10]  21.39  21.38 1.37 1.27  19.11  23.67 1.00     1030     1048
##    y1[1]    6.40   6.34 4.16 3.96  -0.32  13.27 1.01     1206     1380
##    y1[2]    7.97   7.98 4.20 4.08   1.22  14.93 1.00     1369     1959
##    y1[3]    9.67   9.66 4.01 3.76   3.11  16.31 1.00     1531     1689
##    y1[4]   11.43  11.56 3.94 3.84   4.88  17.95 1.00     1285     1835
##    y1[5]   12.95  12.91 3.80 3.60   6.72  19.22 1.00     1803     1812
##    y1[6]   14.85  14.81 3.85 3.74   8.56  21.26 1.00     2093     1891
##    y1[7]   16.42  16.40 3.74 3.62  10.20  22.37 1.00     1891     1798
##    y1[8]   18.01  18.05 3.91 3.75  11.80  24.32 1.00     1861     1922
##    y1[9]   19.76  19.72 3.95 3.78  13.29  26.12 1.00     1891     1991
##    y1[10]  21.37  21.41 3.93 3.72  14.95  27.85 1.00     1925     1672

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -123.73 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -14.5206    0.00175654   0.000469138      0.9852      0.9852       31    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -14.52
##    b0         4.33
##    b1         2.04
##    s          0.67
##    m0[1]     10.44
##    m0[2]     21.25
##    m0[3]      7.28
##    m0[4]     19.43
##    m0[5]     12.53
##    m0[6]     13.09
##    m0[7]     16.77
##    m0[8]     12.53
##    m0[9]     18.81
##    m0[10]     5.59
##    m0[11]     8.52
##    m0[12]     4.62
##    m0[13]    16.78
##    m0[14]    14.95
##    m0[15]    20.58
##    m0[16]    12.81
##    m0[17]    20.58
##    m0[18]    17.46
##    m0[19]    18.76
##    m0[20]    19.97
##    y0[1]     10.53
##    y0[2]     20.86
##    y0[3]      8.25
##    y0[4]     19.23
##    y0[5]     12.43
##    y0[6]     13.62
##    y0[7]     16.10
##    y0[8]     12.95
##    y0[9]     19.39
##    y0[10]     7.30
##    y0[11]     6.65
##    y0[12]     4.40
##    y0[13]    19.75
##    y0[14]    24.10
##    y0[15]    20.91
##    y0[16]    13.48
##    y0[17]    20.55
##    y0[18]    17.68
##    y0[19]    18.38
##    y0[20]    22.28
##    m1[1]      4.62
##    m1[2]      6.47
##    m1[3]      8.32
##    m1[4]     10.16
##    m1[5]     12.01
##    m1[6]     13.86
##    m1[7]     15.71
##    m1[8]     17.55
##    m1[9]     19.40
##    m1[10]    21.25
##    y1[1]    -13.47
##    y1[2]      7.75
##    y1[3]      4.61
##    y1[4]      9.69
##    y1[5]     12.94
##    y1[6]     15.74
##    y1[7]     16.48
##    y1[8]     15.49
##    y1[9]     19.11
##    y1[10]    18.78
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.46 -16.12   1.27 1.03 -19.12 -15.07 1.01      612      686
##    b0       3.96   4.04   0.65 0.64   2.77   4.87 1.01      353      433
##    b1       2.09   2.08   0.11 0.10   1.92   2.28 1.01      409      513
##    s        0.81   0.77   0.24 0.21   0.47   1.23 1.00      752      660
##    m0[1]   10.22  10.26   0.39 0.39   9.53  10.79 1.01      388      557
##    m0[2]   21.31  21.31   0.41 0.38  20.63  21.97 1.00     1096     1025
##    m0[3]    6.98   7.04   0.52 0.51   6.04   7.71 1.01      357      487
##    m0[4]   19.44  19.44   0.34 0.32  18.87  20.01 1.00     1484     1155
##    m0[5]   12.37  12.38   0.32 0.32  11.82  12.86 1.01      461      657
##    m0[6]   12.94  12.95   0.31 0.31  12.42  13.44 1.01      500      664
##    m0[7]   16.72  16.71   0.28 0.27  16.26  17.21 1.00     1466     1157
##    m0[8]   12.36  12.38   0.32 0.32  11.81  12.86 1.01      460      657
##    m0[9]   18.81  18.80   0.33 0.30  18.26  19.35 1.00     1615     1232
##    m0[10]   5.25   5.33   0.60 0.58   4.16   6.08 1.01      354      474
##    m0[11]   8.25   8.31   0.46 0.46   7.42   8.92 1.01      363      528
##    m0[12]   4.26   4.34   0.64 0.62   3.08   5.15 1.01      353      438
##    m0[13]  16.72  16.71   0.28 0.27  16.26  17.21 1.00     1467     1157
##    m0[14]  14.85  14.84   0.28 0.28  14.39  15.33 1.01      792     1264
##    m0[15]  20.62  20.63   0.38 0.35  19.98  21.25 1.00     1225     1005
##    m0[16]  12.65  12.66   0.31 0.31  12.11  13.14 1.01      480      698
##    m0[17]  20.62  20.63   0.38 0.35  19.98  21.25 1.00     1225     1005
##    m0[18]  17.42  17.41   0.29 0.28  16.95  17.93 1.00     1644     1192
##    m0[19]  18.75  18.75   0.32 0.30  18.21  19.30 1.00     1625     1232
##    m0[20]  19.99  19.99   0.36 0.33  19.39  20.57 1.00     1356     1072
##    y0[1]    9.71  10.21  28.05 1.34   3.69  16.25 1.00     1851     1929
##    y0[2]   20.36  21.27  49.97 1.31  15.45  25.83 1.00     1686     1587
##    y0[3]    6.28   6.97  32.37 1.40   2.15  12.18 1.00     1641     2012
##    y0[4]   20.92  19.45  35.11 1.22  14.72  24.54 1.00     1992     1872
##    y0[5]   12.23  12.37  19.04 1.21   7.74  17.49 1.00     1633     1707
##    y0[6]   12.88  12.92  46.89 1.21   8.06  17.92 1.00     1914     1857
##    y0[7]   17.18  16.73  19.70 1.22  12.04  22.59 1.00     2077     1973
##    y0[8]   12.07  12.35  18.45 1.20   7.57  16.81 1.00     1950     1929
##    y0[9]   19.03  18.83  22.33 1.25  13.97  24.27 1.00     2225     1922
##    y0[10]   5.78   5.30  28.75 1.35  -0.74   9.92 1.00     1446     1931
##    y0[11]   8.76   8.25  20.27 1.28   3.08  12.48 1.00     1576     1954
##    y0[12]   3.47   4.27  22.70 1.45  -0.77   9.75 1.00     1574     1963
##    y0[13]  14.83  16.72 102.97 1.20  11.72  20.90 1.00     2058     1960
##    y0[14]  22.15  14.80 278.41 1.18  10.24  20.07 1.00     1562     1763
##    y0[15]  21.28  20.59  33.42 1.30  15.59  25.97 1.00     1871     2024
##    y0[16]  12.31  12.70  31.48 1.19   7.01  17.85 1.00     1933     1828
##    y0[17]  15.48  20.61 258.70 1.36  15.04  26.11 1.00     1948     1859
##    y0[18]  18.33  17.37  48.11 1.30  11.96  22.64 1.00     2147     1974
##    y0[19]  16.65  18.74 108.24 1.20  13.48  23.16 1.00     1967     1861
##    y0[20]  20.10  19.96  26.40 1.28  15.30  26.17 1.00     1922     1927
##    m1[1]    4.26   4.34   0.64 0.62   3.08   5.15 1.01      353      438
##    m1[2]    6.15   6.23   0.55 0.54   5.14   6.92 1.01      355      484
##    m1[3]    8.04   8.11   0.47 0.47   7.20   8.72 1.01      362      507
##    m1[4]    9.94   9.99   0.40 0.40   9.23  10.52 1.01      382      569
##    m1[5]   11.83  11.86   0.34 0.34  11.25  12.35 1.01      436      626
##    m1[6]   13.73  13.73   0.29 0.29  13.25  14.21 1.01      583      977
##    m1[7]   15.62  15.61   0.28 0.27  15.17  16.11 1.01     1033     1280
##    m1[8]   17.52  17.51   0.30 0.28  17.04  18.02 1.00     1659     1122
##    m1[9]   19.41  19.41   0.34 0.32  18.84  19.97 1.00     1489     1127
##    m1[10]  21.31  21.31   0.41 0.38  20.63  21.97 1.00     1096     1025
##    y1[1]    4.88   4.25  15.43 1.48   0.13   9.23 1.00     1223     1960
##    y1[2]    3.27   6.13 121.80 1.37   1.28  11.12 1.00     1255     1943
##    y1[3]    6.97   8.08  62.10 1.35   2.96  13.20 1.00     1551     2031
##    y1[4]   11.30   9.97  40.66 1.32   4.54  15.12 1.00     1539     1915
##    y1[5]   12.60  11.84  19.61 1.20   6.92  16.97 1.01     1280     1595
##    y1[6]   14.87  13.75  57.63 1.22   8.61  18.83 1.00     1771     1979
##    y1[7]   15.24  15.66  32.76 1.28  10.62  21.18 1.00     2010     1788
##    y1[8]   18.06  17.48  43.64 1.18  11.90  22.49 1.00     1890     1931
##    y1[9]   19.62  19.39   9.69 1.24  15.14  24.86 1.00     2063     1973
##    y1[10]  21.49  21.29  22.72 1.24  16.13  26.62 1.00     1828     2004

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.797
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.456
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -6228.25 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       56      -105.906   0.000736808    0.00147673           1           1       66    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -105.91
##    m[1]       4.26
##    m[2]      24.33
##    s[1,1]     5.20
##    s[2,1]    11.48
##    s[1,2]    11.48
##    s[2,2]    47.82
##    xy1[1]     4.13
##    xy1[2]    22.59
##    cr         0.73
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median    sd   mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -101.69 -101.31  1.79  1.61 -105.23 -99.50 1.00      582      771
##    m[1]      4.25    4.23  0.50  0.48    3.48   5.09 1.00      887      861
##    m[2]     24.49   24.69  2.97  2.76   19.34  29.04 1.02      374      445
##    s[1,1]    6.80    6.39  2.17  1.84    4.14  11.10 1.00     1468     1697
##    s[2,1]   14.41   13.52  9.67  8.97    0.77  31.71 1.01      374      479
##    s[1,2]   14.41   13.52  9.67  8.97    0.77  31.71 1.01      374      479
##    s[2,2]   76.64   64.30 44.95 34.81   28.01 163.32 1.01      443      525
##    xy1[1]    4.26    4.22  2.63  2.52   -0.11   8.64 1.00     2016     1955
##    xy1[2]   24.61   24.66  9.17  8.34    9.66  39.38 1.00     1305     1309
##    cr        0.61    0.70  0.27  0.20    0.05   0.89 1.01      406      491

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.607
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -2060.46 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -11.1021   0.000200091   0.000411365       0.962       0.962       38    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -11.10
##    b0        12.87
##    b1         2.04
##    s          2.08
##    m0[1]     23.18
##    m0[2]     20.68
##    m0[3]     21.08
##    m0[4]     28.17
##    m0[5]     18.08
##    m0[6]     27.39
##    m0[7]     18.13
##    m0[8]     18.88
##    m0[9]     24.07
##    y0[1]     23.12
##    y0[2]     21.31
##    y0[3]     21.43
##    y0[4]     27.30
##    y0[5]     23.80
##    y0[6]     28.29
##    y0[7]     19.74
##    y0[8]     17.98
##    y0[9]     28.01
##    m1[1]     12.98
##    m1[2]     14.93
##    m1[3]     16.88
##    m1[4]     18.83
##    m1[5]     20.79
##    m1[6]     22.74
##    m1[7]     24.69
##    m1[8]     26.64
##    m1[9]     28.59
##    m1[10]    30.55
##    y1[1]     12.25
##    y1[2]     13.97
##    y1[3]     16.55
##    y1[4]     18.01
##    y1[5]     21.85
##    y1[6]     19.63
##    y1[7]     28.27
##    y1[8]     25.39
##    y1[9]     30.85
##    y1[10]    29.58
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.15 -11.73 1.50 1.25 -15.07 -10.53 1.00      395      441
##    b0      13.08  12.98 2.79 2.50   8.65  17.94 1.01      258      295
##    b1       2.00   2.02 0.56 0.49   1.07   2.91 1.01      296      380
##    s        2.90   2.65 1.05 0.76   1.74   4.91 1.00      700      880
##    m0[1]   23.23  23.23 1.03 0.91  21.55  24.94 1.00     1539     1106
##    m0[2]   20.78  20.74 1.10 0.98  19.01  22.63 1.00      552      738
##    m0[3]   21.17  21.14 1.06 0.95  19.48  22.94 1.00      633      706
##    m0[4]   28.15  28.15 1.90 1.75  25.12  31.22 1.00      672      878
##    m0[5]   18.22  18.17 1.54 1.32  15.74  20.87 1.01      306      361
##    m0[6]   27.38  27.40 1.72 1.58  24.67  30.14 1.00      754      933
##    m0[7]   18.26  18.21 1.53 1.32  15.80  20.89 1.01      307      361
##    m0[8]   19.00  18.94 1.38 1.17  16.82  21.34 1.00      341      436
##    m0[9]   24.11  24.12 1.11 0.98  22.30  25.95 1.00     1544     1093
##    y0[1]   23.28  23.31 3.11 2.59  18.27  28.11 1.00     1841     1611
##    y0[2]   20.73  20.69 3.31 2.91  15.67  26.05 1.00     1602     1383
##    y0[3]   21.15  21.06 3.19 2.67  16.26  26.41 1.00     1821     1531
##    y0[4]   28.08  28.06 3.72 3.12  21.82  34.24 1.00     1126     1191
##    y0[5]   18.20  18.14 3.50 3.07  12.87  23.90 1.00     1013     1320
##    y0[6]   27.33  27.28 3.40 2.95  21.89  32.96 1.00     1460     1488
##    y0[7]   18.32  18.21 3.52 3.05  12.97  23.81 1.00      914     1216
##    y0[8]   19.02  19.04 3.53 2.95  13.61  24.66 1.00     1028     1373
##    y0[9]   24.12  24.11 3.22 2.83  19.05  29.35 1.00     1940     1736
##    m1[1]   13.19  13.10 2.77 2.49   8.80  18.02 1.01      258      295
##    m1[2]   15.11  15.04 2.27 2.02  11.57  19.07 1.01      264      312
##    m1[3]   17.03  16.98 1.81 1.56  14.25  20.13 1.01      280      368
##    m1[4]   18.96  18.90 1.39 1.19  16.76  21.32 1.00      338      433
##    m1[5]   20.88  20.85 1.09 0.97  19.13  22.68 1.00      571      736
##    m1[6]   22.80  22.80 1.01 0.88  21.18  24.48 1.00     1363     1221
##    m1[7]   24.72  24.72 1.20 1.04  22.79  26.69 1.00     1355     1173
##    m1[8]   26.64  26.65 1.55 1.40  24.16  29.07 1.00      859      992
##    m1[9]   28.57  28.57 2.00 1.86  25.34  31.83 1.00      639      838
##    m1[10]  30.49  30.50 2.48 2.28  26.42  34.56 1.00      543      760
##    y1[1]   13.11  13.13 4.06 3.61   6.44  19.86 1.00      704      956
##    y1[2]   15.13  15.08 3.87 3.37   8.96  21.71 1.00      790     1112
##    y1[3]   16.90  16.81 3.59 3.16  10.98  22.71 1.00      818     1071
##    y1[4]   18.96  19.02 3.39 2.96  13.57  24.40 1.00     1243     1074
##    y1[5]   20.88  20.87 3.28 2.79  15.95  26.06 1.00     1776     1736
##    y1[6]   22.83  22.84 3.23 2.82  17.90  28.00 1.00     1784     1786
##    y1[7]   24.63  24.64 3.25 2.79  19.40  29.66 1.00     1849     1614
##    y1[8]   26.60  26.60 3.37 3.10  21.27  31.82 1.00     1827     1740
##    y1[9]   28.59  28.60 3.71 3.21  22.69  34.35 1.00     1249     1409
##    y1[10]  30.41  30.40 3.89 3.47  24.26  36.72 1.00     1075     1441

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -135.249 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       28      -18.8079   0.000961409   1.54484e-05           1           1       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -18.81
##    b0         8.88
##    b1         3.06
##    s          3.18
##    m0[1]     24.37
##    m0[2]     20.62
##    m0[3]     21.22
##    m0[4]     31.87
##    m0[5]     16.72
##    m0[6]     30.70
##    m0[7]     16.79
##    m0[8]     17.92
##    m0[9]     25.70
##    y0[1]     25.42
##    y0[2]     21.83
##    y0[3]     20.08
##    y0[4]     33.79
##    y0[5]     17.84
##    y0[6]     30.79
##    y0[7]     13.15
##    y0[8]     18.34
##    y0[9]     23.62
##    m1[1]      9.05
##    m1[2]     11.98
##    m1[3]     14.91
##    m1[4]     17.85
##    m1[5]     20.78
##    m1[6]     23.71
##    m1[7]     26.64
##    m1[8]     29.57
##    m1[9]     32.50
##    m1[10]    35.43
##    y1[1]      8.34
##    y1[2]      8.50
##    y1[3]     17.36
##    y1[4]     20.74
##    y1[5]     16.89
##    y1[6]     25.68
##    y1[7]     25.23
##    y1[8]     29.82
##    y1[9]     36.11
##    y1[10]    37.43
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -19.50 -19.12 1.50 1.26 -22.58 -17.82 1.01      392      416
##    b0       6.96   7.52 3.79 3.23   0.72  11.74 1.02      237      258
##    b1       3.50   3.37 0.78 0.62   2.54   4.86 1.01      280      353
##    s        4.49   4.11 1.58 1.21   2.69   7.56 1.01      566      682
##    m0[1]   24.68  24.56 1.41 1.25  22.61  27.22 1.00     1137     1045
##    m0[2]   20.40  20.46 1.45 1.25  18.05  22.64 1.00      427      637
##    m0[3]   21.08  21.14 1.40 1.24  18.84  23.29 1.00      492      695
##    m0[4]   33.27  32.79 2.68 2.28  29.86  38.48 1.01      575      484
##    m0[5]   15.93  16.13 2.05 1.75  12.51  18.77 1.01      268      360
##    m0[6]   31.93  31.50 2.43 2.10  28.81  36.61 1.01      640      508
##    m0[7]   16.01  16.20 2.03 1.74  12.60  18.84 1.01      269      360
##    m0[8]   17.30  17.44 1.83 1.56  14.24  19.95 1.01      290      385
##    m0[9]   26.21  26.04 1.54 1.35  24.03  29.04 1.00     1207      871
##    y0[1]   24.70  24.80 5.04 4.28  16.59  32.66 1.00     2019     1682
##    y0[2]   20.33  20.46 5.07 4.49  12.08  28.26 1.00     1811     1642
##    y0[3]   21.00  21.05 4.90 4.25  12.87  28.89 1.00     2100     1716
##    y0[4]   33.31  32.92 5.34 4.80  25.42  42.14 1.00     1340     1389
##    y0[5]   15.78  16.02 5.25 4.70   6.96  23.68 1.00     1249     1659
##    y0[6]   31.96  31.49 5.46 4.80  24.08  41.18 1.00     1438     1354
##    y0[7]   15.95  16.16 5.26 4.36   7.55  24.04 1.00     1024     1131
##    y0[8]   17.20  17.40 4.92 4.53   8.81  24.70 1.00     1049     1289
##    y0[9]   26.23  26.05 4.96 4.42  18.78  34.61 1.00     1789     1600
##    m1[1]    7.16   7.72 3.75 3.20   0.98  11.89 1.02      237      252
##    m1[2]   10.51  10.92 3.06 2.64   5.52  14.52 1.01      239      271
##    m1[3]   13.87  14.14 2.41 2.06   9.85  17.13 1.01      249      284
##    m1[4]   17.22  17.37 1.84 1.57  14.13  19.87 1.01      288      384
##    m1[5]   20.58  20.65 1.44 1.24  18.28  22.82 1.00      442      644
##    m1[6]   23.93  23.84 1.37 1.21  21.89  26.38 1.00     1014     1034
##    m1[7]   27.28  27.05 1.67 1.48  24.95  30.47 1.01     1120      786
##    m1[8]   30.64  30.23 2.20 1.90  27.79  34.89 1.01      729      558
##    m1[9]   33.99  33.50 2.82 2.37  30.40  39.44 1.01      546      497
##    m1[10]  37.35  36.75 3.50 2.89  32.90  44.06 1.01      459      447
##    y1[1]    7.20   7.69 5.91 5.29  -2.96  15.64 1.00      550      841
##    y1[2]   10.48  11.02 5.58 4.82   1.15  18.57 1.00      711      951
##    y1[3]   13.99  14.27 5.42 4.77   5.01  21.89 1.00      859     1073
##    y1[4]   17.22  17.51 5.03 4.44   9.02  24.86 1.00     1498     1452
##    y1[5]   20.73  20.65 4.85 4.52  12.88  28.49 1.00     2050     1752
##    y1[6]   23.83  23.80 4.80 4.13  16.50  31.64 1.00     1760     1554
##    y1[7]   27.10  26.97 4.93 4.33  19.82  35.34 1.00     1638     1523
##    y1[8]   30.68  30.37 5.27 4.48  22.47  39.55 1.00     1452     1151
##    y1[9]   33.90  33.40 5.47 4.72  26.16  43.74 1.00     1328      891
##    y1[10]  37.38  36.79 5.98 5.28  28.51  47.85 1.00      866      815

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -26.386 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       -11.803     0.0014423   8.29677e-05           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -11.80
##      p        0.81
##      se       0.23
##      sp       0.43
##      ppv      0.64
##      npv      0.11
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -17.90 -17.59 1.29 1.10 -20.38 -16.44 1.00      731      895
##      p      0.51   0.52 0.29 0.38   0.06   0.95 1.00      684      598
##      se     0.27   0.25 0.11 0.12   0.10   0.47 1.00     2710     1570
##      sp     0.44   0.44 0.15 0.16   0.19   0.70 1.00     2426     1502
##      ppv    0.39   0.32 0.29 0.34   0.02   0.91 1.00      722      548
##      npv    0.40   0.35 0.29 0.34   0.02   0.91 1.00      720      585

2x2 cross table

Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p01)
n11~B(n1.,p11)

ex16-1.stan

data {
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p01;
  real<lower=0,upper=1> p11;
}
model {
  n01~binomial(n0,p01);
  n11~binomial(n1,p11);
}
generated quantities {
  real RR;
  RR=p11/p01;
  real OR;
  OR=(p11/(1-p11))/(p01/(1-p01));
}
n0=20
n01=rbinom(1,n0,0.3)
n1=20
n11=rbinom(1,n1,0.6)
data=list(n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -36.372 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       -22.217    0.00761187   0.000346553           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -22.22
##      p01      0.15
##      p11      0.55
##      RR       3.67
##      OR       6.92
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -26.62 -26.28 1.08 0.70 -28.81 -25.64 1.01      744      998
##      p01    0.18   0.17 0.08 0.08   0.07   0.34 1.00     1182      823
##      p11    0.54   0.55 0.10 0.10   0.37   0.71 1.01     1856     1223
##      RR     3.83   3.17 2.89 1.56   1.45   8.36 1.01     1198      975
##      OR     8.03   5.97 8.60 4.02   1.79  20.13 1.01     1315     1157

bimodal distribution, mixed normal distribution

n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -247 100  6 -521 -524
## 
## Clustering table:
##  1  2  3 
## 35 31 34
rst$parameters
## $pro
## [1] 0.346 0.308 0.346
## 
## $mean
##      1      2      3 
## -5.104 -0.269  5.189 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.953
## 
## 
## $Vinv
## NULL
plot(rst)

ex17-1.stan

data {
  int N;
  vector[N] y;;
}
parameters {
  simplex[3] h; //ratio of mix
  real m1;
  real m2;
  real m3;
  real<lower=0> s1;
  real<lower=0> s2;
  real<lower=0> s3;
}
model {
  s1~cauchy(0,5);
  s2~cauchy(0,5);
  s3~cauchy(0,5);
  for (i in 1:N) {
    vector[3] p;
    p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
    p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
    p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
    target+=log_sum_exp(p);
  }
}
mdl=cmdstan_model('./ex17-1.stan')

data=list(N=n,y=y)
mle=mdl$optimize(data=data)
## Initial log joint probability = -707.624 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       39      -246.361   0.000153551    0.00108076           1           1       56    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -246.36
##      h[1]     0.34
##      h[2]     0.34
##      h[3]     0.32
##      m1       5.23
##      m2      -5.13
##      m3      -0.26
##      s1       0.92
##      s2       0.87
##      s3       1.13
mcmc=goMCMC(mdl,data,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 2 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 3 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 4 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 3 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 2 finished in 1.7 seconds.
## Chain 3 finished in 1.7 seconds.
## Chain 4 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 4 finished in 1.7 seconds.
## Chain 1 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 1 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 1 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 1 finished in 7.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3.0 seconds.
## Total execution time: 7.1 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median    sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -255.32 -253.80  6.00 2.17 -275.94 -251.20 1.05       56       11
##      h[1]    0.36    0.34  0.09 0.05    0.27    0.58 1.05       48       11
##      h[2]    0.33    0.33  0.05 0.05    0.25    0.41 1.03      111     2246
##      h[3]    0.31    0.33  0.09 0.05    0.04    0.41 1.09       29       11
##      m1      2.47    5.12  4.45 0.35   -5.27    5.49 1.57        6       27
##      m2     -0.13   -0.29  3.67 4.16   -5.29    5.36 2.41        4       25
##      m3    -18.26   -4.92 78.70 1.08  -17.04    0.01 1.81        8       11
##      s1      1.10    0.97  0.55 0.16    0.76    2.66 1.08       31       11
##      s2      1.11    1.05  0.30 0.23    0.77    1.64 1.25       11      115
##      s3      2.01    1.02 12.54 0.24    0.76    2.02 1.29       10       17